Cracking Condition of Cohesionless Porous Materials in Drying Processes 
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The invasion of air into porous systems in drying processes is often localized in soft materials, 
such as colloidal suspensions and granular pastes, and it typically develops in the form of cracks 
before ordinary drying begins. To investigate such processes, we construct an invasion percolation 
model on a deformable lattice for cohesionless elastic systems, and with this model we derive the 
condition under which cracking occurs. A Grifhth-like condition characterized by a dimensionless 
parameter is proposed, and its validity is checked numerically. This condition indicates that the 
ease with which cracking occurs increases as the particles composing the material become smaller, 
as the rigidity of the system increases, and as the degree of heterogeneity characterizing the drying 
processes decreases. 

PACS numbers: 46.50.+a,64.70.fm,81.05.Rm,83.80.Hj 
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I. INTRODUCTION 

Contraction that results from drying often causes the 
formation of cracks in wet granular materials and pastes 
of colloid suspensions. Such materials are soft due to 
weak cohesion among constituent particles, and they pos- 
sess unique properties characterizing the crack formation 
that they exhibit, including memory effects, slow crack 
growth and a diversity of crack patterns [ll-lllj. Consid- 
erable effort has been devoted to investigating the for- 
mation of crack patterns in the contexts of physics, soil 
mechanics and geology, while cracking conditions have 
been investigated mainly in engineering applications. 

During a drying process, the first crack in a paste gen- 
erally appears in a capillary state, in which all pores (i.e. 
the spaces between solid particles) are filled with liquid, 
and the cracking process consists essentially of the in- 
vasion of air into these pores. We believe that cracking 
conditions can be deduced directly from material proper- 
ties related to drying, because on the microscopic level, 
cracking is indistinguishable from ordinary drying in the 
elementary processes of air invasion. However, drying 
and cracking have been treated separately in previous 
theories, although some studies have investigated the ef- 
fect of large deformation and the order of ordinary drying 
and cracking in drying processes [l^ - [l6j . 

The main goal of this paper is to determine the Griffith 
criterion associated with the drying process of cohesion- 
less porous systems. For this purpose, we focus on the 
limit of slow drying in a uniform elastic system. The slow 
drying limit can be realized by decreasing the relative hu- 
midity h quasistatically during the drying processes with 
fixed temperature T and atmospheric pressure P. Under 
such conditions, both the water distribution and elastic 
deformation of the material can be assumed to be in ther- 
mal equilibrium or quasi-equilibrium. We also assume 



that the elastic relaxation is much faster than the redis- 
tribution of liquid, and we ignore plasticity, although it 
is believed to be important for many type of pastes. 

To investigate the systems of interest, we extend an in- 
vasion percolation (IP) model to include elastic interac- 
tions. It has been established that IP models are faithful 
models of drying processes in porous materials, despite 
the great simplification they employ of treatin g liquid 
distributions as binary distributions on a lattice |17H25j . 
Our extended model describes the crack-like invasion of 
air in soft systems with large rigidity. With it, we de- 
termine the condition for the formation of the first crack 
in terms of the free energy of the system. This condi- 
tion is determined from the elastic properties of the ma- 
terial, the heterogeneity characterizing the drying pro- 
cesses, and the size of particles forming the paste. 

In this paper, we investigated a two-dimensional sys- 
tem corresponding to a cross-section of a uniform layer 
of paste. We regard the bottom surface to be a fixed 
boundary, and the top surface to be a free boundary, 
and assume that the evaporation of liquid occurs only 
from the top surface. This paper is constructed as fol- 
lows. We propose a theoretical model based on thermo- 
dynamic considerations in Sec. |TT1 In Sec. IIIIl we report 
the results of numerical simulations using this model. We 
study the conditions governing the invasion of air and the 
formation of the first crack in Sec. II VI and fVl respectively. 
Conclusions are given in Sec. IVII 
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We regard the system as a mixture of solid particles 
and liquid forming a paste layer, and the environment as 
air containing vapor, which exists both inside and outside 
the paste layer. We treat this system and environment 
thermodynamically. In order to find the thermal equilib- 
rium state for a given h with fixed (T, P), we introduce a 
free energy J = F + P{Vs + Vi)-Hy (T, P, h)Ni , where F, 
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Vs, fj-v, Ni and Vi are the Helmholtz free energy, the total 
volume of the constituent particles, the chemical poten- 
tial of the vapor, the number of liquid molecules, and the 
volume of the liquid, respectively. We ignore the effect 
of gravity and assume that Vs and Vi/Ni are constant. 

The Helmholtz free energy, F, is given by the sum 
F = Fe + Fi + Fi. Here, F^ represents the deforma- 
tion energy of the system, which results from interactions 
among particles. The surface energy, Fi, is determined 
by the liquid distribution and increases as the invasion 
of air progresses, as shown in Appendix \X\ The free en- 
ergy of the liquid, Fi, is given hy Fi = -PiVi + fiiNi, 
as obtained from the Gibbs-Diihem relation, where Pi is 
the hydrostatic pressure of the liquid, and the chemical 
potential of the liquid, /i/ , is identical to ■ Substituting 
these forms into the above expression for J, and defining 
p = P — Pi, we obtain 



J = Fe + Fi+ p{T, P, h)Vi + const. 



(1) 



The equilibrium state of the paste for given h is deter- 
mined by minimizing J. 

The difference between the hydrostatic pressure of the 
fluid in the pores and the atmospheric pressure, p, is 
called the negative pore pressure in soil mechanics. The 
pressure p is uniform throughout the system in equilib- 
rium states, and determined uniquely by h for given T 
and P, in accordance with the Kelvin equation. We 
adopt p as the control parameter in place of h, because 
p corresponds to the driving force of air invasion. In the 
processes we consider, p increases in time, while h de- 
creases as a function of p. 

Here we give a brief remark concerning the minimiza- 
tion of J for processes occurring at fixed p. First, we 
note that, even when p increases quasistatically, a region 
into which air has invaded often expands abruptly, ex- 
hibiting behavior similar to that in an avalanche process. 
However, if the redistribution and desiccation of liquid 
are sufhciently slow, such a process can still be regarded 
as proceeding slowly, and therefore we can assume that 
the system is approximately in mechanical equilibrium 
throughout. While Vi decreases gradually, instantaneous 
equilibrium states (characterized by instantaneous values 
of Vi) can be obtained by minimizing F. Technically, this 
implies that we can always use the minimum principle of 
J by including Vi into the set of control parameters dur- 
ing such a process, because the difference between J and 
F, P{Vs + Vi)- /^^,(r, P, h{T, P,p))Ni, is constant if Vi is 
fixed, in addition to {T,P,p). 



A. Rigid materials 

We use a regular lattice composed of M cells each of 
volume AV to represent the system. The liquid distri- 
bution on the lattice is described by the set of variables 
{</>!, •• • , (/)m}, where (j)m represents the "dryness" of the 
mth cell. Each cell is assumed to be either wet {(j)m = 0) 
or dry (0^ = 1), as in percolation models. In order for 



such a treatment to be valid, we must assume that the 
cells possess microscopic volumes. We adopt the length 
of a cell as the unit of length in our model, which is 
proportional to the linear extent of a particle, r, fixing 
all other microscopic properties. In addition, we assume 
that all wet cells have the same liquid volume fraction, 
fiu, and that no liquid exists in dry cells, for simplic- 
ity. The volume occupied by solid particles in each cell 
is (1 - Vyj)AV. 

The total volume of liquid and surface energy are 



M 



and 



F,, 



m=l 



M 

Jla AAr, 



const. 



(2) 



(3) 



respectively, where 7/a is the surface tension of the liquid- 
air interface. The quantity jiaAAm is the increase in the 
surface energy needed to change the state of the mth cell 
from wet to dry. This represents the resistance to drying. 
In order to investigate heterogeneous material, we assume 
that A Am varies among the cells, with an average value 
AA. The characteristic pressure in a drying process, 



P7 = lla 



AA Jia 



VwAV r 



(4) 



is typically on the order of the surface tension divided by 
r. 

Because Fg is constant for rigid systems, we have 

M 

J = ^ [yiaAAmCpm. + PVwi'^ - 4'7n)AV] + coust. (5) 
m— 1 

Introducing the quantities J = J/v^ip^, 7™ = AAm/ AA 
and p = p/p-y, this expression can be written in the sim- 
plified form 



J 



M 

E 

m— 1 



AV{jrn4'm + p(l " 0m)}, 



(6) 



after removing the constant from J. 

This free energy is minimized when the liquid distri- 
bution satisfies the conditions 



(t),n = for 7„i > p 

(f>jn = 1 for 7m < p 



(7) 



If these conditions are applied only to cells adjacent to 
the dry- wet interface, this process corresponds to that de- 
scribed by the conventional IP model without a trapping 
rule jl7Hl9| . In the present work, we assume that 7m is 
distributed uniformly over the interval [1 — Aj, 1 -I- Aj] 
and that Aj ^ 1. For large systems, air invades when 
p ~ 1 — A^ from the drying surface, and it percolates 
above a certain threshold p ^ Pc < 1 + Aj, as is well 
known. 
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B. Extension to elastic materials 

We regard each cell as an elastic tile subject to uniform 
strains, in order to describe the deformation of the sys- 
tem. The elastic energy Fe is determined by the strains 



as F,{U„.} = E™=i ^VfeiUm), where f/„ 



- t (™) 

"'■aP 



)is 



the strain tensor of the mth cell and fe is the free energy 
density. In order to simplify the situation, we make the 
following assumptions: (al) all cells have the same elas- 
tic properties; (a2) fe does not depend on 0^; (a3) AA^ 
does not depend on Um ■ Heterogeneity is introduced into 
the system only through the drying process. 

Coupling of (f>„i and U„i is introduced through Vi . The 
dilation of a cell results from the influx of liquid for a wet 
cell and an influx of air for a dry cell. To account for such 
phenomena, we add the volume expansion ratio ui™"* to 
Vw in Eq. This yields the free energy 

M 



J = Fe{U,n} + J2^1laAA, 



m—1 

+p{vw + - 'l)m)AV] + const, (8) 

where Einstein's summation rule is applied to repeated 
Greek suffixes. Then, introducing Um = Um/vw and 

fe{Um) = feiUjn)/VwP'y, WC obtain 
M 

J = ^Z\y{/e(t/„O-f7m0m+P(l+4a^)(l-0m)}- (9) 
m—1 

We note that the assumption (aS) can be weakened 
slightly in the case that AA^ depends on Um linearly 
as A Am = AA{'-/m + where 7' is a constant, be- 

cause an expression identical to Eq. (|9]) can be obtained 
through an appropriate transformation of the variables. 

With the elastic energy included, the conditions to de- 
termine the liquid distribution are revised from those ap- 
pearing in Eq. ^ to 



= 1 for 7„i <pil+ Uaa 



= for 7m > p{l + wi"^, 

- ~{m)\ 



(10) 



These conditions imply that the resistance to drying de- 
creases as a cell expands, because the expansion of a cell 
decreases the cost in surface energy required to remove a 
unit volume of liquid. Here, we note that the change of 
a state from wet to dry is reversible, and may be caused 
simply by the redistribution of liquid, rather than desic- 
cation, as observed by L. Xu et.al [26.]- Also, desiccation 
in the vicinity of the free surface generally induces liquid 
flow which causes the invasion of air far from this surface. 

The free energy has a minimum with respect to {Um} 
in the equilibrium state. Minimizing Eq. ([9]) in the con- 
tinuum limit gives the stress balance equation 



da, 



af) 



dxj 



0, 



(11) 



where Xj represents the space coordinates and the 
stresses dap are defined by 



dfe 

du. 



ap 



'P{1 - (t>)Sap. 



(12) 



Wet regions are subject to compressive pressure p from 
both the free surface and the interface with dry regions. 



C. Elastic energy of cohesionless materials 

Soft materials that exhibit drying cracks, such as col- 
loid suspensions and wet granular materials, generally 
have nonlinear elastic properties. In many cases, the co- 
hesive interactions of the constituent particles are very 
weak in comparison with excluded volume interactions. 
The elastic energy and moduli practically vanish unless 
the system is subject to compressive stresses. Therefore, 
we assume that fe{U) = for Uaa > and that the 
elastic moduli vanish for Uaa — 0. 

We need to choose an appropriate function of fe{U) 
for Uaa < 0, because there is no established general 
constituent relation. If we assume an isotropic ana- 
lytic function for the elastic moduli, fe can be approx- 
imated in the form of a third-order elasticity as fe = 
—^w^a ~ fo'^ Small Uaa, where A and are 

positive constants. Goehring reported that this third- 
order elasticity accurately describes the results of com- 
pression tests with cornstarch paste [13] ■ Another choice 

is fe = \/—u^i{^Uaa + f^^ap)' which is obtained theo- 
retically by assuming Hertzian contacts among particles 
and affine deformation [l^ . 

We assume the following general homogeneous form 
including these choices: 



fe{U) = giu^i) 



^Kula + G (uap - ^U^rjSaP 



(13) 

where K and G are positive constants. The function g{x) 
takes the power-law form 



aix) 



{-x)"-'^ x<0 
x>0 ' 



(14) 



where 1/^2 for third-order elasticity and = 3/2 in the 
case of Hertzian contacts. The bulk modulus and rigidity 
depend on the state of the material, due to the nonlin- 
earity. As shown in Appendix |b1 they are proportional 
to g{uaa)K and g{uaa)G, respectively, for isotropic com- 
pressive states, and 2G < v{y + \)K is required for ther- 
mal stability. We investigated our model mainly for large 
as Eq. ([T^ generally holds for small deformations. 



III. NUMERICAL SIMULATIONS 

We consider a uniform layer from whose top surface 
(which is a one-dimensional interface) liquid is desic- 
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FIG. 1. Triangular lattice with Nx x Ny vertices and M — 
2Ny{Nx — 1) cells. The a;-axis is normal to the surfaces of the 
system, with the value of x representing the distance from the 
top surface. Dry cells are enclosed by bold lines. 
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FIG. 2. The drying process following the crack-like invasion 
of air. The graph displays the results for the dry fraction (p 
and the maximum depth of a dry cell, Xmax, as functions 
of p obtained from a numerical simulation with {K, G) — 
(10^2 X 10*), A-y = 0.001 and ^ Ny ^ 80. The bottom 
snapshots correspond to (p — 1)/ A'y = 1.19 and 2.79. 



cated. We carried out numerical simulations using third- 
order elasticity (j/ = 2) on a two-dimensional triangular 
lattice, as shown in Fig. [1] The lattice has x Ny 
vertices, and the layer thickness is H = \/3(Nx — l)/2. 

In order to avoid erroneous numerical convergence due 
to the singularity at Uaa = for g{x), we used the smooth 
function 







for 6a; < —1 
for \bx\ < 1 
for 6a; > 1 



(15) 



with a large positive constant b 
form given in Eq. (jl4[) . 



10 , instead of the 



The numerical method we used is essentially the same 
as that used in [2^. The deformation of the lat- 
tice is described by the displacements of the vertices, 
{iti, • • • , itAf} {N = NxNy), which determine {Um}. The 
top surface is a free boundary in contact with air. The in- 
terface between air and wet cells is treated as the dry- wet 
interface. The bottom surface is a fixed boundary with 
respect to {un} and a reflecting boundary with respect 
to {(t>m}- Periodic boundary conditions are used along 
the y direction. The following procedures were repeated 
in the numerical simulations from the initial conditions 
in which all cells were wet and undeformed. 

(PI): The displacements {m„} were calculated by mini- 
mizing J in Eqs. (151) and (|13p using the conjugate 
gradient method [30|. 

(P2): The conditions (|T0| were checked for all cells con- 
tacting the dry-wet interface. If these conditions 
were satisfied, we increased p by Ap and then re- 
turned to (PI). 

(P3): If the mth cell did not satisfy these conditions, 
we changed 0m from wet (dry) to dry (wet) and 
returned to (PI). 

In (P3), if we found more than one cell that did not sat- 
isfy the conditions , we changed only the state of the 
most unstable cell, i.e., that with the largest deviation 
from the condition 7^ = p{l + u^cIm)- Under this proce- 
dure, Vi changes gradually for fixed p. In the simulations 
whose results are presented here, we used Ap = Z\7/500, 
and the tolerance of the conjugate gradient method was 
3 X 10-". 

Our numerical results indicate that the crack-like in- 
vasion of air occurs readily in soft systems with large 
rigidity and small heterogeneity. 

Figure [2] depicts a typical process of crack-like inva- 
sion and subsequent drying. As soon as a cell on the 
top surface dries, air penetrates rapidly into the bulk 
and this results in the formation of a one-dimensional 
dry region. This dry region expands gradually from 
both the top surface and the crack line as p increases. 
The graph in this figure plots the fraction of the entire 
system in the dry region, (f> = X]m=i 't'm/M, and the 
maximum depth of a dry cell divided by the thickness, 
Xmax = maxm{a;m|0m = ^} / H , as functions of p. It is 
found that soft material resists the invasion of air through 
shrinkage. When shrinkage occurs, the pressure required 
for air invasion to begin is larger than that in the case 
of the conventional IP model, p ~ 1 — Aj. On the other 
hand, the percolation threshold decreases drastically if 
the crack-like invasion of air occurs. 

Figure [3] displays typical snapshots for 9 sets of {K, G) 
and fixed A^. These were taken after a dry region had 
developed beyond a depth of approximately half of the 
total depth. For sufficiently large K or small G, the 
air invasion process occurs in the same manner as for the 
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FIG. 3. Snapshots of air invasion for A-y = 0.001, A" = 10^ 10*, 10^ and G/K = 0.02, 0.1, 0.5. In all cases, = Ny = 80. The 
number appearing in each figure indicates the value of (p ~ 1)/ A'y for that snapshot. The deformation of the lattice is scaled 
by a factor of 5 for viewability, and the volume expansion ratio Uaa is indicated by the gray scale, which corresponds to the 

interval [0^ - AU, Ui + AU], where Ui is given by Eq. ((20} and AU = 1/(6^10^). 



conventional IP model. It alvifays occurs at G = for any 
value of K, as described below. The dry region expands 
intermittently and gradually as the pressure p increases. 
Contrastingly, for large G and small K, crack-like air 
invasion occurs first. This dry region essentially corre- 
sponds to a mode-I crack, because the cells in this region 
are expanded (i.e. Uaa > 0), while the surrounding wet 
region shrinks. For sufficiently large G and small K, 
however, shear bands (niode-II cracks) often form ahead 
of cracks in the wet region. The directions of the shear 
bands reflect the anisotropy of the triangular lattice. The 
wet cells in such shear bands expand until Uaa — 0, and 
crack-like air invasion develops along shear bands. Shear 
bands sometimes appear and disappear during the inva- 
sion of air. 

Figure S] elucidates the dependence of the heterogene- 
ity for fixed {K,G). As increases, the air invasion 
process changes from crack-like to conventional IP-like. 
The value of Aj at which this change occurs depends on 
{K,G), as discussed in Sec. IVl 



IV. INVASION OF AIR 

When there is no crack, the system is uniform and 
contracts only along the x direction as p increases with 
drying. Such uniaxial compressive states are described 
by (f) — and 

M=(^^^^ l)^Ua, (16) 

where Ua is determined from the free boundary condition 
i^xx — dfe/dUa +J5 = 0. Substituting these values into 
Eq. (fTi)) . we obtain 

Ua^- I — r and feitfa) = —pUa- 

(17) 

If the top surface is sufficiently long, the invasion of air 
begins when the second condition in (fTO| holds at a cell 
on the surface for which we have 

I- A-f^p{l + Ua), (18) 
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FIG. 4. Snapshots of air invasions for Ay = 0.001, 0.01, 0.1 and fixed {K, G) = (lO", 10^). Tiiese snapshots are displayed in the 



same manner as those in Fig. [S] In all cases, A^^; 
the respective values of Ay. 



Ny = 80. The number in each figure indicates the value of (p — I)/ A-y for 



because the smallest jm is 1 — Ay. The pressure at which 
the invasion of air begins is determined by this condition 
and increases as K+G decreases. Because the right-hand 
side of Eq. ([T8|) has a maximum at Ua = —1^/(1^+ 1), the 
invasion of air never occurs for systems so soft that the 
condition 

K + G> 2(^1 + ^^ (1 - Zi7) ~ 2 (^1 + (19) 

holds. This case corresponds to wet sintering in this sys- 
tem 0,liM5i. 

Let us first consider a drying process for G — 0, that 
is, the case with no rigidity. In this case, the system 
does not exhibit cracking and dries in the same manner 
as in the conventional IP model. We can solve the stress 
balance equation (ITT|) easily in the case G = 0, and we 
obtain isotropic stress states with fr^^ = and 

for all wet cells and Uaa > for dry cells. The elastic 
energy of a wet cell is the same as that of a cell in the 
isotropic compressive state: 

fem = '^PU^, where t/, ^ ^ J J ) . (21) 

The pressure at which the mth cell is allowed to dry 
is determined from y„i — p(l -I- Ui{p)), in accordance 
with the conditions in (jlOp. along with the condition 
Ui > —vjiy -f 1). Because the pressure always increases 
monotonically as 7„i increases, the order of drying does 
not depend on K. 

The above conclusion concerning the nature of the dry- 
ing can also be obtained by considering the free energy, 
J. To show this, we first note that additional energy is re- 
quired to dry cells on an elastic lattice, because the elastic 



energy stored in a cell is dissipated quickly through di- 
lation following drying. In order to develop a dry region 
D for a fixed value of p, the amount by which the free 
energy J decreases, must be larger than the dissipation, 

i?^ ^ Zil//e(f/(fy)), (22) 

where TJm^^'^ is the strain on the mth cell at the time 
that it changes from wet to dry. This condition can be 
written 

J(p; <i>) - J(p; D)-R=Y,AV (p(l + C/,) - 7™) > 0, 

(23) 

because fe{Um'^^^) — fe{Ui), and Eq. ^ can be rewrit- 
ten as J{p-D) = T.,ntD^y{hm+P{l + U,)] + 

Sms-D^^'T'"' where D = (p corresponds to the initial 
state, in which all cells are wet. When the invasion of 
air begins with the condition (|18p . the right-hand side is 
negative, because p(l -I- IJi) = 1 — Ay < y^- Therefore, 
the dry region cannot expand quickly, and hence there is 
no cracking. 

V. CRACKING 

Next, we consider the case G 7^ and investigate the 
cracking condition. Cracking in our system consists of a 
quasi-one-dimensional invasion of air. It is facilitated by 
the release of energy from a surrounding wet region. Just 
as in the situation considered in the previous section, a 
crack develops if the amount by which the free energy 
J decreases is larger than the dissipation R when the 
invasion of air begins. 

If there appears a crack of length L that is perpen- 
dicular to the free surface in a sufficiently large system 
with fixed p, this crack forms an air-invaded region D 
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and it causes the elastic energy stored in D to dissipate 
locahy. The wet region S surrounding the crack, from 
which the elastic energy is released, has an area of ap- 
proximately L X L, because L is the only characteristic 
length in this system. The states in this region change 
from the uniaxial compressive state, characterized by Ua, 
to an approximately isotropic compressive state, charac- 
terized by Ui. Measuring L by the number of cells in D, 
the decrease of J is estimated as 

JiP;^)- J{p;D) 

~ L'^AV {feiUa) - +p{Ua - tj{)] 

+ LAv[uUa)+p{l + Ua)-lD{L)Y (24) 
where the average of 7„i in D is represented by 
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1D{ 



(25) 



meD 



The elastic energy of the mth cell in D, fe{Um), dis- 
sipates with drying at — p(l + u^aa)- A cell is de- 
formed at a crack tip and expands by at least an amount 

wL™^ — Ua — (7m — 1 -|- A^)/p = 0{A^) in comparison 
with Eq. ([T5)) . However, as shown in Appendix [Cl the 
energy in this case is the same as fe{Ua) to first order in 
Aj, because the cell contacts a dry-wet interface. Thus, 
the minimum dissipation for air invasion is 



R = LAvUUa)+0{iA-ff). 



(26) 



The cracking condition J{p; 4 
p and D is approximated as 



J(p;L>)-i?> Ofor fixed 



L 



{feiUa)- fem+PiUa-Ui)} 
>lD{L)~p{l + Ua)^lD{L) 



1 + Z\7, (27) 



where we have used Eq. (|T8| . The right-hand side of this 
expression is positive and order A^^. The cracking condi- 
tion of the system is determined by finding the minimum 
value of ^d{L) over all possible crack paths D. This 
yields 

L [hiila) - fem+p{Ua " U.)] > c{L) A^ , (28) 

where c(L)Z\7 = min^j (7/3 (L) — 1 + Aj). The quantity 
c{L) increases from c(l) ~ to a value less than 1 as L 
increases. 

The condition derived above corresponds to the Grif- 
fith criterion for the first crack in a drying process. It 
indicates that there is a critical crack length L = Lc be- 
yond which unstable crack growth occurs. The left-hand 
side of Eq. (^5)1 is the energy released when the crack 
advances by one cell. The right-hand side, c{L)A'-f, is 
the additional energy required to dry a cell at a crack 
tip. The critical length, Lc, is determined by the non- 
dimensional parameter 



r ; 



^7 {feiUa) - feiUi)+p{Ua - U,)}' 



(29) 
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FIG. 5. Dependence of AX on Xo for K = 1000 and A-y = 
0.001. The results for each value of G were calculated from 
the average of Xi over 8 numerical simulations on a lattice 



with Nx 



Ny =40. The error bars for the plots with G = 20 



and 200 are omitted for clarity. 



as Lc = c{Lc)T. Substituting Eqs. ^ and ^ 

into Eq. we obtain 



1 



A-1 



K-A-f 



p{Ua-U.) g.{§)p'+i 



(30) 



where 

9u{x) 



v+l) 2 



[l-{l + x)-^]. (31) 



The quantity gu{x) is an increasing function of x and 
is approximately proportional to a;/2 for a; <C 1. From 
Eq. ([HI), we find p ~ 1 for > 1 and A-^ < 1. 

In order to verify the criterion appearing in Eq. (j28p 
on the basis of our numerical simulations, it is necessary 
to quantify the transition from the invasion of the con- 
ventional IP model to crack-like invasion. We calculated 
the quantities 



1 



M 

E 



X^(f>n 



(32) 



for k 
AX = 



— and 1, and investigated the dependence of 
^2X1 - X§ on Xo. 
The conventional IP model is characterized by the self- 
affine growth of a dry- wet interface. Assuming the in- 
terface to be a single-valued function x = h{y) in the 
initial stages, the standard deviation a/ {{h — {h)y) in- 
creases as a power function of the average height (h), 
where (•) represents the average over y. Because 4>„i — 1 
for < Xm < h(jj), the standard deviation corresponds to 
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FIG. 6. Dependence of a' on 1/G for various K and 
A"/. These results were obtained by applying a least-square 
method to a long-log plot of AX and Xo for Xo < 1 to cal- 
culate a'. 



FIG. 8. L/c(L) for a triangular lattice of = Ny = 40. The 
quantity c{L) used here is the average over 100 sets of {7m} 
created from different random seeds. The inset displays an 
example of a one-dimensional path with length L = 11. 
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FIG. 7. Dependence of a' on F. The data in Fig. |6] are 
replotted with respect to F, defined in Eq. (|30p . 



AX, as we have Xi ~ {^^ dxx) = {h'^)/2 and Xq ~ (h). 
Contrastingly, in a crack-like process, \/Xi and, hence, 
AX are proportional to Xq, because in such a process, 
a one-dimensional dry region develops. Figure [5] displays 
a typical dependence of AX on Xq. As seen there, the 
invasion changes from that described by the conventional 
IP model to crack-like invasion as G increases. The ex- 
ponent a', defined by AX cx Xq at ~ 0, is approxi- 
mately 0.5 for the invasion of the conventional IP model 
and 1 for crack-like invasion. 

Figure H] displays the exponent a' as a function of 1 /G 
for various K and Aj. The change of a' from 1.0 to 0.5 
corresponds to the transition of the invasion type. The 



data are replotted with respect to F in Fig. [71 It is seen 
that all data approximately collapse to a single curve for 
large K, and the transition occurs at approximately the 
same value of F in each case, near F ~ 10. This result is 
consistent with the criterion appearing in Eq. (j28p and 
indicates that Lc/c{Lc) — 10 in our simulations. 

The function c(L) is determined from the distribution 
function of 7™ and the lattice properties. Investigating 
all one-dimensional paths of length L, extending from the 
top surface in the x direction, we determined c{L) from 
the minimum value of the sum of 7™ over a path. Fig- 



ure [8] displays L/c{L) for the Nx 



40 triangular 



lattice we used. It is seen that, because c{L) increases 
with L from ~ at L = 1, L/c{L) first decreases and 
then increases as L increases from 1. For the portion 
of this graph in which L/c{L) is increasing, the value of 
L along the curve is L^., the crack length beyond which 
a crack grows unstably. For the portion of this graph 
in which LjciV) is decreasing, the value of L along the 
curve corresponds to the depth of the dry region that de- 
velops following the first invasion due to heterogeneity. 
The minimum of c{L)/L corresponds to the transition 
point of F below which Lc vanishes and the first inva- 
sion spontaneously induces an initial crack that develops 
into crack-like invasion. The value found here for F is 
consistent with F ~ 10, obtained from Fig. [71 

The quantity F generally depends on the particle size. 
In the original, unsealed system, {K, G) correspond to 
{K,G) = {p-y/v'^){K,G). For fixed elastic properties, 
represented by {K, G) and heterogeneity, represented by 
zi7, we find from Eq. ^ that F given in Eq. ([HOjl with 
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1 increases with the size of a particle as 



k) 



9u 



(33) 



This implies that cracking occurs only when the particle 
size is small. Similar results have been reported previ- 
ously in some experiments |26i] . 

We now discuss the validity of the cracking criterion 
given in Eq. p8)) . This criterion seems valid for large 
r, even when shear bands appear in advance of the first 
crack. This is because this crack will simply follow the 
path of the shear band. Most of the difference in the 
free energy is consumed when the shear band forms, be- 
cause a shear band enables compression in the surround- 
ing region, due to the large strain it creates. By contrast, 
cracks tend to become wide for sufficiently soft materials. 
This is due to the fact that if a one-dimensional crack 
were to develop in such materials, the large expansion 
at a crack tip would cause additional dissipation, larger 
than 0((Z\7)2), in Eq. The dependence of a on T 

seen in Fig. [7] deviates from the master curve for small 
if, because we assumed a crack to be a one-dimensional 
region when we derived our criterion. In order to gener- 
alize our criterion to be applicable to such cases, we need 
to evaluate elastic strains at the tip of a fat crack. 

We have investigated the first crack appearing in a 
uniform layer without initial cracks or flaws, except mi- 
croscopic heterogeneity of drying properties. If the layer 
contains a initial crack or a macroscopic air-invaded re- 
gion initially, it can develop and divide the system at 
smaller values of p. When p is very small, the last term 
in Eq. (|24)) can be approximated as LAV{—jd{L)}, and 
the dissipation term can be ignored. In this case, the 
cracking condition can be written 

L {/e(f/a) - fc{U,)+p{Ua - U^)} > Jd{^), (34) 

because ^d{L) is approximately constant for i ^ 1. As- 
suming L to be the layer thickness, iJ, this condition 
gives 



7d(oo) 



9> 



H 



(35) 



as the smallest value of the pressure for which cracking 
can occur. For the original, unsealed system, we have 
p oc H~'^ , and p does not depend on r, because H is 
scaled by the unit of length, which is proportional to r. 
These dependences are consistent with the experimental 
results obtained by Man et.al. [l6j . 



of paste-like materials. We derived a cracking condition 
that applies to cohesionless porous systems taking the 
same form as the Griffith criterion, after eliminating lo- 
cal dissipation accompanied by air invasion. The Griffith 
energy corresponds to an additional energy required for 
drying, not the surface energy of the liquid-air interface 
itself. We find that crack-like air invasion occurs for soft 
materials with larger rigidity and less heterogeneity in 
the properties characterizing the drying process. Also, 
this criterion explains why cracking does not occur for 
systems composed of large particles. 

For systems in which there is fast drying or large de- 
formation, the cracking condition will difi^er from that 
derived here, because in such situations, there are compli- 
cations that were not accounted for in the present work. 
Specifically, in the case of fast drying, the pore pressure 
will become non-uniform, while in the case of larger de- 
formation, plastic deformation will appear. While it is 
important to elucidate such phenomena, these problems 
are beyond the scope of the present work. 
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Appendix A: Equilibrium conditions 

Let us consider thermal equilibrium states of a paste 
system for given (T,P,h), where the relative humidity, 
h, is given hy h = Py/ P* , when Py is the vapor pressure, 
and P* is its saturated value. We assume that both the 
particles and the liquid composing the system are incom- 
pressible and that the vapor is an ideal gas, for simplicity. 

A condition of mechanical equilibrium, Laplace's law, 
states that the liquid pressure. Pi, and the atmospheric 
pressure, P, are related as 



P^ P - Pi = JlaK. 



(Al) 



Here k is the mean curvature, which is constant every- 
where on the liquid-air interface in an equilibrium state. 

The negative pore pressure p is determined from h by 
the Kelvin condition, 



VI. CONCLUSIONS 

We proposed an invasion percolation model for a co- 
hesionless elastic material to investigate drying processes 



p = -ksTpi \ogh, 



(A2) 



for an ideal gas, where fc^ is Boltzmann's constant and 
pi = Ni/Vi is the number density of liquid molecules. 
This equation is derived from the chemical equilibrium 
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condition fii (T, Pi) = Hv (T, Py ) and the condition for a 
fiat interface, ^i(T,P) — fj.y{T,P*), which provides the 
definition of P* . In the derivation of Eq. (|A2|) , we have 
used the relations dpii /dPi — 1/pi and fj,y = fc^T log P^ + 
const. 

Another mechanical equilibrium condition, Young- 
Dupre's law, 

Jla COS 6* = - ^sl , (A3) 

holds at the contact points of the liquid-air interface 
and the surfaces of solid particles. Here, 9 is the con- 
tact angle. The interface energy of the paste is the 
sum of the surface energies of the liquid-air, solid-air, 
and solid-liquid interfaces. Explicitly, we have Fi = 
jiaAia + jsaAsa + "fsiAsi, where Amn and jmn are the 
surface area and the surface tension of the interface indi- 
cated by their suffixes, respectively. As the total area of 
the solid surfaces, Asa + Asi, is approximately constant, 
Fi can be written as 

= 7ia^ + const., (A4) 

where A = Aia + cos 9 Asa, after substituting Eq. (jA3p . 
The invasion of air causes A and thus Fi to increase. 



Appendix B: Isotropic compressive states 

The free energy J of a wet cell is minimal with respect 
to Uai3 in the isotropic compressive state. For Eqs. ()13|) 
and l|14p . J depends on as 

(Bl) 

When the strain tensor deviates from that of an isotropic 
state by [/q^, taking the form — Ui -|- Uap), 

this equation can be approximated to second order in 
Uap as 

fe P^aa 

~ const. + |p - !^±^{-U.,rK^ U,U^^ 



Applying the condition that this quantity be minimal 
gives Eq. ^ and {v + l)vk > 2G. The bulk modulus 
and rigidity in the isotropic compressive state are {v + 
\)v{-U^Y-^k/2 and {-tj.y-^G, respectively. 



Appendix C: Perturbation from a uniaxial 
compressive state 

A wet cell in contact with a dry-wet interface has the 
same stress conditions, Uxx = <^xy = 0, as the unixaial 
compressive state, described by Ua given in Eq. (|16|) . 
where the x-axis is perpendicular to the interface, while 
the stress (jyy depends on the volume expansion ratio. 

The elastic energy is /e(C/) = fe{Ua)+(7eap{Ua)5uap + 
0{5U'^) for U = Ua + SU, where aeap = dfe/duafs = 
o'ap — pSai3 ■ The first-order term can be rewritten as 



— (^"^xx^^exx ~^ ^xy^^exy ^" ^yy^^eyy) i,^^) 



because fe = o'eajSUap / {v + 1) for the homogeneous 
function given in Eq. (|14l) and dfe = d'ea/sdua/S = 
[deafiduap + ddeapUap) / {1 + v) . Equation ()Cip van- 
ishes, because Sd^xx = Sdexy — for fixed p, due to 
the stress conditions and Uyy{Ua) = 0. Thus, we have 
MU) - feiUa) ^ 0{SU^) and 

(TeapSilaiS = -pSUxx + <^eyySUyy = 0. (C2) 

The quantity 5U is determined from Eq. (|C2|) and Sa^y = 
for a given Suaa = Suxx + ^Uyy, and thus we find 

/e(t/) - UUa) = 0{{SUa.c.f). (C3) 
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